home *** CD-ROM | disk | FTP | other *** search
/ Mac Mania 4 / MacMania 4.toast / / Demo's / Igor Demo Pro / 3 PutContentsIn Igor Pro Folder / Technical Notes / Igor Tech Notes / TN020-B Peak Areas / Peak Areas Folder / procedure < prev    next >
Text File  |  1993-09-24  |  25KB  |  1,079 lines
  1. | Peak Areas Experiment V1.3 - 9/24/93
  2. Menu "Macros"
  3.     m_so
  4.     Submenu m_b
  5.         m_bi
  6.         m_ba
  7.         m_bd
  8.         m_bf
  9.         m_s
  10.         m_bs
  11.     End
  12.     "-"
  13.     Submenu "Identify And Measure Peaks in Steps"
  14.         m_ii
  15.         "Auto Identify Peaks…/3"
  16.         "Ignore Peaks Between Cursors/4"
  17.         "Identify Maxima At Cursor A/5"
  18.         "Identify Minima At Cursor B/6" 
  19.         "Construct Baseline From Peaks.../7"
  20.         "-"
  21.         "Measure Identified Peaks..."
  22.     End
  23.     "Auto Identify And Measure Peaks..."
  24.     "Zoom To Peak.../8"
  25.     "-"
  26.     Submenu m_a
  27.         m_ai
  28.         m_aa
  29.         m_ala
  30.     End
  31.     "-"
  32.     m_ah
  33.     m_cug
  34.     m_cw
  35. End
  36.  
  37.  
  38. Macro InitializeMostEverything(w,wx,wb)
  39.     String w=g_w,wx=g_wx,wb=g_b
  40.     Prompt w,p_w,popup, WaveList("*",";","")
  41.     Prompt wx,p_wx,popup,calc+WaveList("*X", ";","")
  42.     Prompt wb,p_b,popup, none+WaveList("*base*",";","")
  43.  
  44.     Silent 1;PauseUpdate
  45.     g_w=w;g_wx=wx; SBs(wb)
  46.     Redimension/D $w
  47.     SameLen(w,breg);$breg=NaN;SameLen(w,areg);$areg=NaN
  48.     Mk(ampsY,2);Mk(ctrsX,2);Mk(widsX,2);Mk(ctrsP,2);Mk(edgsP,4);
  49.     Mk(areaNB,2);Mk(areaX1,2);Mk(areaX2,2);
  50.     Mk(bdcw,2);Mk(pdcw,2)
  51.     if(exists(wx)==1)
  52.         Redimension/D$wx
  53.     endif
  54.     if(exists(wb)==1)
  55.         Redimension/D $wb
  56.     endif
  57.     AppWv(w,wx,"");HideAreaTags()
  58. End
  59.  
  60. Macro InitBaseLineFit(w,wx)
  61.     String w=g_w,wx=g_wx
  62.     Prompt w,p_w,popup, WaveList("*",";","")
  63.     Prompt wx,p_wx,popup,calc+WaveList("*X", ";","")
  64.     
  65.     Silent 1;PauseUpdate
  66.     g_w=w;g_wx=wx
  67.     SameLen(w,breg);AppWv(breg,wx,"");Modify mode($breg)=1;$breg=NaN;AppCrs(w)
  68. End
  69. Macro AddRegionToFit()
  70.     Silent 1;PauseUpdate
  71.     CheckTwoCursors(g_w);$breg[V_start,V_theEnd]=$g_w[p]
  72. End
  73. Macro DeleteRegionFromFit()
  74.     Silent 1;PauseUpdate
  75.     CheckTwoCursors(g_w);$breg[V_start,V_theEnd]=NaN
  76. End
  77.  
  78. Macro FitBaselineAtRegions(w,wx,wr,fit)
  79.     String w=g_w,wx=g_wx,wr=fbar_wr,fit=fbar_fit
  80.     Prompt w,p_w,popup, WaveList("*",";","")
  81.     Prompt wx,p_wx,popup,calc+WaveList("*X", ";","")
  82.     Prompt wr,"region wave (baseline weighting wave)",popup, none+breg
  83.     Prompt fit,"baseline function",popup,S_funcs
  84.     
  85.     Silent 1;PauseUpdate
  86.     g_w=w;g_wx=wx;fbar_wr=wr;fbar_fit=fit
  87.     ChkLen(w,wx);ChkLen(w,wr)
  88.     String ow="W_BaselineFit",wtmp="W_tmp",opts="";SameLen(w,ow)
  89.     Variable t=0,typ=floor(strsearch(S_funcs,fit,0)/10)
  90.     
  91.     if(exists(wx)==1)
  92.         opts+="/X="+wx
  93.     endif
  94.     if(exists(wr)==1)
  95.         Dup(wr,wtmp)
  96.         $wtmp=(numtype($wtmp)==0)
  97.         opts+="/W="+wtmp
  98.         ar_ex=3
  99.     else
  100.         ar_ex=1
  101.     endif
  102.     t=str2num(fit[7,8])-1
  103.     String command="CurveFit "+fit[1,7]+", $w"+opts
  104.     Execute command
  105.     KillWv(wtmp)
  106.     Mk(bcw,2);$bcw={NaN,K0,K1,K2,K3,K4};Redimension/N=(2+t) $bcw
  107.     Mk(bdcw,2);$bdcw={1,0,numpnts($w)-1,typ,t};Dup(bdcw,"W_PM")
  108.     if(exists(wx)==1)
  109.         $ow=PolyMorph($bcw,$wx[p])
  110.     else
  111.         $ow=PolyMorph($bcw,x)
  112.     endif
  113.     AppWv(ow,wx,"")
  114.     SBs(ow)
  115.     ar_wfit=ow
  116. End
  117.  
  118. Macro SubtractBaseline(w,wx,wb,ow)
  119.     String w=g_w,wx=g_wx,wb=rb_b,ow=rb_ow
  120.     Prompt w,p_w,popup, WaveList("*",";","")
  121.     Prompt wx,p_wx,popup,calc+WaveList("*X", ";","")
  122.     Prompt wb,p_b,popup, "_From Fit Peaks_;"+WaveList("*base*",";","")
  123.     Prompt ow,"output baseline-corrected wave",popup,new+WaveList("*",";","")
  124.  
  125.     Silent 1;PauseUpdate
  126.     g_w=w;g_wx=wx
  127.     ChkLen(w,wx);ChkLen(w,wb)
  128.     String cw,dcw,tcw="W_tmp0",tdcw="W_PM"
  129.     if(exists(ow)!=1)
  130.         NewWv(w,"NoBase");ow=S_Wave
  131.     else
  132.         SameLen(w,ow)
  133.     endif
  134.     rb_ow=ow
  135.     if(exists(wb)!=1)
  136.         cw=bcw
  137.         dcw=bdcw
  138.         if(numtype($dcw[1])!=0)
  139.             Abort "Run Fit Peaks first!"
  140.         endif
  141.         Variable terms=$dcw[4],s=$dcw[1],e=$dcw[2],funcs=$dcw[0]
  142.         Dup(cw,tcw);Redimension/N=(2+terms) $cw
  143.         Dup(dcw,tdcw);$tdcw[0]=1
  144.         wb="W_tmp";SameLen(w,wb);$wb=NaN
  145.         if(exists(wx)==1)
  146.             $wb[s,e]=PolyMorph($tcw,$wx[p])
  147.         else
  148.             $wb[s,e]=PolyMorph($tcw,x)
  149.         endif
  150.     endif
  151.  
  152.     $ow=$w-$wb
  153.     AppWv(ow,wx,"");Modify zero(left)=1
  154.     KillWv(wb);KillWv(tcw)
  155.     g_w=ow
  156.     SBs("_None_")
  157. End
  158.  
  159. Macro InitAreaBetweenCursors(what,w,wx,wb,anno)
  160.     Variable what=mir_what,anno=mir_anno
  161.     String w=g_w,wx=g_wx,wb=g_b
  162.     Prompt what,"integration type",popup, mir_pwh
  163.     Prompt w,p_wd,popup, WaveList("*",";","")
  164.     Prompt wx,p_wx,popup,calc+WaveList("*X", ";","")
  165.     Prompt wb,p_b,popup,none+WaveList("*base*", ";","")
  166.     Prompt anno,"area annotation",popup,none+"_one tag on graph_;_one tag per peak_"
  167.     
  168.     Silent 1;PauseUpdate
  169.     mir_what=what;g_w=w;g_wx=wx;SBs(wb);mir_anno=anno
  170.     ChkLen(w,wx);ChkLen(w,wb);AppWv(w,wx,"");AppWv(wb,wx,"")
  171.     SameLen(w,areg);AppWv(areg,wx,"");Modify mode($areg)=1;$areg=NaN;AppCrs(w)
  172.     HideAreaTags();Mk(areaNB,2);Mk(areaX1,2);Mk(areaX2,2)
  173. End
  174. Macro AreaBetweenCursors()
  175.     Silent 1;PauseUpdate
  176.     String w=g_w,wx=g_wx,wb=g_b,text,wv=w
  177.     CheckTwoCursors(w)
  178.     Variable/D s=V_start,en=V_theEnd,wi=0,oi,wlx,wrx
  179.     if (mir_anno!=3)
  180.         $areg=NaN;HideAreaTags()
  181.     endif
  182.     $areg[s,en]=$w[p]
  183.     iterate (2)
  184.  
  185.         wlx=pnt2x($wv,s);wrx=pnt2x($wv,en)
  186.         oi=wi
  187.         AreaKind(mir_what,wv,wx,s,en)    | outputs V_Area,S_Text
  188.         wi=V_Area
  189.         Print wv+" "+S_Text+" = "+num2str(wi)
  190.         if(exists(wb)!=1)
  191.             break
  192.         endif
  193.         wv=wb
  194.     loop
  195.     text=S_Text
  196.     if(exists(wb)==1)
  197.         Print w+" - "+wb+" "+text+" = "+num2str(oi-wi)
  198.         wi=oi-wi;text += "\r(less baseline)"
  199.     endif
  200.     if(mir_anno!=1)
  201.         Tag/C/N=$("area"+num2istr(g_atgn))/A=MC $w,(pnt2x($w,s)+pnt2x($w,en))/2, text+"\r= "+num2str(wi)
  202.         if(mir_anno==3)
  203.             g_atgn+=1
  204.         endif
  205.     endif
  206.     InsertPoints 0,1,$areaNB,$areaX1,$areaX2;$areaNB[0]=wi;$areaX1[0]=V_startX;$areaX2[0]=V_theEndX
  207.     Sort $areaX1,$areaX1,$areaX2,$areaNB
  208. End
  209.  
  210. Macro HideAreaTags()
  211.     do
  212.         Tag/K/N=$("area"+num2istr(g_atgn))
  213.         g_atgn-=1
  214.     while (g_atgn>=0)
  215.     g_atgn=0
  216. End
  217.  
  218. Macro ListAreas()
  219.     DoWindow/F AreaReportTable
  220.     if(V_Flag==0)
  221.         AreaReportTable()
  222.     endif
  223. End
  224.  
  225. Macro AutoIdentifyAndMeasurePeaks()
  226.     InitIdentifyPeaks()
  227.     ShowOnlyDataAndBase()
  228.     AutoIdentifyPeaks(2,g_w,g_wx,g_b,g_bx,1)
  229.     String kn
  230.     if(tfp_pol==1) | positive peaks
  231.         kn=wmnp | baseline from neg pks
  232.     else
  233.         kn=wmxp
  234.     endif
  235.     ConstructBaselineFromPeaks(g_w,g_wx,kn)
  236.     MeasureIdentifiedPeaks(mir_what,pfr_title,pfr_sort,pfr_order,tfp_pol)
  237. End
  238.  
  239. Macro InitIdentifyPeaks(w,wx,wb,pol,box,what)
  240.     String w=g_w,wx=g_wx,wb=g_b
  241.     Variable box=g_bx,pol=mpr_pol,what=mir_what
  242.     Prompt w,p_w,popup, WaveList("*",";","")
  243.     Prompt wx,p_wx,popup,calc+WaveList("*X", ";","")
  244.     Prompt wb,p_b,popup, none+WaveList("*base*",";","")
  245.     Prompt pol,"peak polarity", popup,"positive;negative"
  246.     Prompt box,p_bx
  247.     Prompt what,"default integration type",popup, mir_pwh
  248.  
  249.     
  250.     Silent 1;PauseUpdate
  251.     g_w=w;g_wx=wx;SBs(wb)
  252.     box=abs(box);mpr_pol=pol;tfp_pol=pol;g_bx=box;mir_what=what
  253.  
  254.     ChkLen(w,wx);ChkLen(w,wb)
  255.     SameLen(w,wmnp);SameLen(w,wmxp);$wmnp=NaN;$wmxp=NaN
  256.     Mk(ampsY,2);Mk(ctrsX,2);Mk(ctrsP,2);Mk(widsX,2);Mk(edgsP,4)
  257.     DoWindow/F PeakFitGraph
  258.     if(V_Flag==0)
  259.         AppWv(w,wx,"");DoWindow/C PeakFitGraph
  260.     else
  261.         AppWv(w,wx,"");HideAreaTags()
  262.     endif
  263.     AppWv(wb,wx,"")
  264.     AppCrs(w);SetAxis/A left;SetAxis/A bottom
  265. End
  266. Macro IdentifyMaximaAtCursorA()
  267.     if(strlen(CsrWave(A))==0)
  268.         Abort "expecting cursor A in target graph on "+g_w+" or similar wave"
  269.     endif
  270.     if(exists(g_wx)==1)
  271.         FindPX(g_w,g_wx,hcsr(A))
  272.         $wmxp[V_peakP]=$g_w[V_peakP]
  273.     else
  274.         $wmxp(hcsr(A))=$g_w(hcsr(A))
  275.     endif
  276. End
  277. Macro IdentifyMinimaAtCursorB()
  278.     if(strlen(CsrWave(B))==0)
  279.         Abort "expecting cursor B in target graph on "+g_w+" or similar wave"
  280.     endif
  281.     if(exists(g_wx)==1)
  282.         FindPX(g_w,g_wx,hcsr(B))
  283.         $wmnp[V_peakP]=$g_w[V_peakP]
  284.     else
  285.         $wmnp(hcsr(B))=$g_w(hcsr(B))
  286.     endif
  287. End
  288. Macro IgnorePeaksBetweenCursors()
  289.     Silent 1;PauseUpdate
  290.     String w=g_w,wx=g_wx,wb=g_b
  291.     CheckTwoCursors(g_w)
  292.     $wmnp[V_start,V_theEnd]=NaN
  293.     $wmxp[V_start,V_theEnd]=NaN
  294. End
  295.  
  296. | Find min,max peaks by analyzing peak data
  297. Macro AutoIdentifyPeaks(what,w,wx,wb,box,extent)
  298.     String w=g_w,wx=g_wx,wb=g_b
  299.     Variable what=tfp_what,box=g_bx,extent=tfp_extent
  300.     Prompt what,"do what with found peaks?",popup, p_aipw
  301.     Prompt w,p_w,popup,WaveList("*", ";","")
  302.     Prompt wx,p_wx,popup,calc+WaveList("*X", ";","")
  303.     Prompt wb,p_b,popup,none+WaveList("*base*", ";","")
  304.     Prompt box,p_bx
  305.     Prompt extent,"extent of Y wave to search",popup,"entire wave;between cursors"
  306.  
  307.     Silent 1;PauseUpdate
  308.     tfp_what=what;g_w=w;g_wx=wx;SBs(wb)
  309.     box=abs(box)
  310.     g_bx=box;tfp_extent=extent
  311.  
  312.     ChkLen(w,wx);ChkLen(w,wb)
  313.     Variable s=0,en=numpnts($w)-1
  314.     if(extent==2) | use cursors
  315.         CheckTwoCursors(w);s=V_start;en=V_theEnd
  316.     endif
  317.     XBoxToPBox(box,w,wx,s,en);box=g_pbox;Print "minima/maxima window is ",box," points..."
  318.     String tw="W_tmp"
  319.     if(exists(wb)==1)
  320.         SameLen(w,tw);$tw[0,s]=NaN;$tw[en,numpnts($tw)-1]=NaN
  321.         $tw[s,en]=$w[p]-$wb[p]
  322.         w=tw
  323.     endif
  324.     SameLen(w,wmnp);SameLen(w,wmxp);
  325.     if(what==2)
  326.         $wmnp=NaN;$wmxp=NaN
  327.     endif
  328.     | Find min peaks, max peaks
  329.     $wmnp[s,en]=IsBestMinMax($wmnp,$w,box/2,p,0)    | 1 if local min or max, else 0
  330.     $wmxp[s,en]=IsBestMinMax($wmxp,$w,box/2,p,1)
  331.     $wmnp[s,en]=$w[p] * sqrt(-1+2*$wmnp[p])    | wave val or NaN
  332.     $wmxp[s,en]=$w[p] * sqrt(-1+2*$wmxp[p])
  333.     AppWv(wmnp,wx,"");AppWv(wmxp,wx,"")
  334.     Modify mode($wmnp)=8,marker($wmnp)=1,mode($wmxp)=3,marker($wmxp)=8
  335.     WaveStats/Q/R=[s,en] $wmxp;Print v_npnts," maxima found"
  336.     WaveStats/Q/R=[s,en] $wmnp;Print v_npnts," minima found"
  337.     KillWv(tw)
  338. End
  339.  
  340. Function IsBestMinMax(wmm,w,win,pt,doMax)
  341.     wave/D wmm,w
  342.     Variable win,pt,doMax
  343.     Variable isMM=IsLocalMinMax(w,win,pt,doMax)
  344.     if(isMM %& (pt > 0))
  345.         Variable p1=limit(pt-win,0,pt-1)
  346.         if(mean(wmm,pnt2x(wmm,p1),pnt2x(wmm,pt-1))>0)    | another min or max already found is too close
  347.             isMM=0
  348.         endif
  349.     endif
  350.     return isMM
  351. End
  352.  
  353. |returns truth that point pt in wave w is local maxima or minima over +/- win points
  354. Function IsLocalMinMax(w,win,pt,doMax)
  355.     wave/D w
  356.     Variable win,pt,doMax
  357.     Variable isMM=1,p=max(0,pt-win),lim=min(numpnts(w)-1,pt+win)
  358.     Variable/D mn=w[pt],mx=w[pt]
  359.     do
  360.         if(doMax)
  361.             if(w[p]>mx)
  362.                 isMM=0
  363.                 break;
  364.             endif
  365.         else
  366.             if(w[p]<mn)
  367.                 isMM=0
  368.                 break;
  369.             endif
  370.         endif
  371.         p+=1
  372.     while(p<=lim)
  373.     return isMM
  374. End
  375. Macro ConstructBaselineFromPeaks(w,wx,knots)
  376.     String w=g_w,wx=g_wx
  377.     String knots=cbfp_pol
  378.     Prompt w,p_w,popup,WaveList("*", ";","")
  379.     Prompt wx,p_wx,popup,calc+WaveList("*X", ";","")
  380.     Prompt knots,"peaks to generate baseline from", popup,wmxp+";"+wmnp
  381.  
  382.     Silent 1;PauseUpdate
  383.     String wb="W_BaselineFit" | output wave
  384.     cbfp_pol=knots;g_w=w;g_wx=wx;SBs(wb)
  385.     ChkLen(w,wx);Dup(w,wb);ChkLen(w,knots)
  386.     WaveStats/Q $knots
  387.     Variable npks,nots=V_npnts
  388.     if(nots<2)
  389.         Abort "need 2 or more knots in knots wave "+knots
  390.     else
  391.         Print "using",nots,"knots"
  392.     endif
  393.     |Create straight line baseline using peaks as knots
  394.     if(exists(wx)==1)
  395.         npks= KnotsXYToBase($w,$wx,$wb,$knots)
  396.     else
  397.         npks= KnotsYToBase($w,$wb,$knots)
  398.     endif
  399.     AppWv(wb,wx,"")
  400. End
  401. | inputs are W_MinKnots and W_MaxKnots
  402. | outputs are
  403. |
  404. | W_PkCenters[i]            Final peak center value after baseline correction
  405. | W_PkAmps[i]                Final peak amplitude after baseline correction
  406. | W_PkFWHMs[i]            Full Width at Half Maximums of the baseline-corrected peaks
  407. | W_PkAreas[i]                Areas of the baseline-corrected peaks
  408. | W_PkX1[i]                Starting x coordinate of area computation. NEW
  409. | W_PkX2[i]                Ending x coordinate of area computation.    NEW
  410. | W_PkXwL                    Starting x coordinate of fwhm (add width to get end x) NEW
  411.  
  412. Macro MeasureIdentifiedPeaks(what,title,srt,order,pol)
  413.     Variable what=mir_what
  414.     String title=pfr_title, srt=pfr_sort,
  415.     Variable order=pfr_order,pol=tfp_pol
  416.     Prompt what,"integration type",popup, mir_pwh
  417.     Prompt title,"report title"
  418.     Prompt srt,"sort report by",popup,"center;amplitude;width;area"
  419.     Prompt order,"sort order",popup,"ascending;descending"
  420.     Prompt pol,"peak polarity", popup,"positive;negative" 
  421.         
  422.     Silent 1;PauseUpdate
  423.     pfr_title=title;pfr_sort=srt;pfr_order=order;tfp_pol=pol;mir_what=what
  424.     String center="W_PkCenters",amplitude="W_PkAmps",width="W_PkFWHMs",area="W_PkAreas"
  425.     String wx1="W_PkX1",wx2="W_PkX2",fc="W_PkXwL",pkw,bsw
  426.     if(pol==1)
  427.         pkw=wmxp
  428.         bsw=wmnp
  429.     else
  430.         pkw=wmnp
  431.         bsw=wmxp
  432.     endif
  433.     WaveStats/Q $pkw
  434.     Variable npks=V_npnts,np1=npks-1,nrows=max(2,npks)
  435.     String text=num2istr(npks)+" peaks"
  436.     Mk(center,nrows);Mk(amplitude,nrows);Mk(width,nrows);Mk(area,nrows);
  437.     Mk(wx1,nrows);Mk(wx2,nrows);Mk(fc,nrows)
  438.     if(npks<1)
  439.         Abort "no peaks"
  440.     else
  441.         print npks," peaks to measure..."
  442.     endif
  443.     Variable pc=numpnts($pkw),p1,p2,lim=numpnts($bsw)-1
  444.     Variable/D m,yb,ab,xd,y0,f1,f2,totalArea=0
  445. Variable ti=ticks
  446.     iterate (npks)
  447.         pc=PrevNonNaN($pkw,pc) | center of right-most peak
  448.         p1=limit(PrevNonNan($bsw,pc),0,lim)    | start
  449.         p2=limit(NextNonNan($bsw,pc),0,lim)    |end
  450.         y0=$g_w[pc]
  451.         if(exists(g_wx)==1)
  452.             $center[i]= $g_wx[pc]
  453.             $wx1[i]= $g_wx[p1]
  454.             $wx2[i]= $g_wx[p2]
  455.         else
  456.             $center[i]= pnt2x($g_w,pc)
  457.             $wx1[i]= pnt2x($g_w,p1)
  458.             $wx2[i]= pnt2x($g_w,p2)
  459.         endif
  460.         xd=$wx2[i]-$wx1[i]
  461.         if(p1==p2)
  462.             m=0
  463.         else
  464.             m=($g_w[p2]-$g_w[p1])/xd
  465.         endif
  466.         yb=$g_w[p1]+m*($center[i]-$wx1[i])
  467.         $amplitude[i] = y0-yb
  468.          AreaKind(what,g_w,g_wx,p1,p2)    | outputs V_Area,S_Text
  469.         if(what==1) |count
  470.             ab=(p2-p1+1)*($g_w[p1]+$g_w[p2])/2
  471.         else
  472.             ab=xd*($g_w[p1]+$g_w[p2])/2
  473.         endif
  474.         V_Area-=ab
  475.         $area[i]= V_Area
  476.         totalArea+=V_Area
  477.         FindLevelXY(g_w,g_wx,pc,p1,(y0+$g_w[p1])/2);f1=V_LevelX
  478.         FindLevelXY(g_w,g_wx,pc,p2,(y0+$g_w[p2])/2);f2=V_LevelX
  479.         $fc[i]=f1
  480.         $width[i]=f2-f1
  481.         if(mod(i+1,25)==0)
  482.             print npks-i,"..."
  483.         else
  484.             printf "%d...",npks-i
  485.         endif
  486. |print i,yb,ab,f1,f2
  487.     loop
  488.     Printf "\rTotal Area=%8.8g\r",totalArea
  489. ti=(ticks-ti)/60
  490. print ti/npks,"seconds per peak"
  491.     DoWindow/F PeakReportTable
  492.     if(V_Flag==0)
  493.         PeakReportTable()
  494.     endif
  495.     text+="\r("+S_Text+")"
  496.     String cmd="Sort "
  497.     if(order==2)
  498.         cmd+="/R  "
  499.     endif
  500.     cmd+="$"+srt+",$center,$amplitude,$width,$area,$wx1,$wx2,$fc"
  501.     Execute cmd        | Sort...
  502.     
  503.     DoWindow/F PeakReportLayout
  504.     if(V_Flag==0)
  505.         PeakReportLayout()
  506.     endif
  507.     ReplaceText/N=info text
  508.     ReplaceText/N=title "\JC"+title
  509.     Modify rows(PeakReportTable)=nrows
  510. End
  511.  
  512.  
  513. Macro ZoomToPeak(peak,anno)
  514.     Variable peak=ztp_pn,anno=mir_anno
  515.     Prompt anno,"area annotation",popup,none+"_one tag on graph_;_one tag per peak_"
  516.  
  517.     Silent 1;PauseUpdate
  518.     ztp_pn=peak;mir_anno=anno
  519.     WaveStats/Q W_PkCenters
  520.     Variable/D npks=V_npnts,x1,x2,xc
  521.     DoWindow/F PeakFitGraph
  522.     if( (V_Flag==1) %& (peak==limit(peak,0,npks-1)) )
  523.         String ex="W_ShowEstWidthsX",ey="W_ShowEstWidthsY",text
  524.         x1=W_PkXwL[peak];x2=x1+W_PkFWHMs[peak]
  525.         Mk(ex,4);Mk(ey,4);$ex={x1,x1,x2,x2}
  526.         x1=W_PkX1[peak];x2=W_PkX2[peak];xc=W_PkCenters[peak];
  527.         SetAxis bottom,x1,x2
  528.         if(exists(g_wx)==1)
  529.             FindPX(g_w,g_wx,x1);x1=pnt2x($g_w,V_PeakP)
  530.             FindPX(g_w,g_wx,x2);x2=pnt2x($g_w,V_PeakP)
  531.             FindPX(g_w,g_wx,xc);xc=pnt2x($g_w,V_PeakP)
  532.         endif
  533.         $ey={0,$g_w(x1)+$g_w(xc),$g_w(x2)+$g_w(xc),0};$ey/=2
  534.         WaveStats/Q/R=(x1,x2) $g_w
  535.         SetAxis left V_min,V_max
  536.         AppWv(ey,ex,"");Modify mode($ey)=1
  537.         Cursor/P A,$ey,1;Cursor/P B,$ey,2
  538.         if(mir_anno!=3)
  539.             HideAreaTags()
  540.         endif
  541.         if(mir_anno!=1)
  542.             sprintf text,"\JCpeak %d\rarea=%g\rfwhm=%g",peak,W_PkAreas[peak],W_PkFWHMs[peak]
  543.             Tag/C/N=$("area"+num2istr(g_atgn))/A=MC $g_w,xc, text
  544.             if(mir_anno==3)
  545.                 g_atgn+=1
  546.             endif
  547.         endif
  548.     else
  549.         Beep;SetAxis/A left;SetAxis/A bottom
  550.     endif
  551. End
  552.  
  553.  
  554. Proc FindLevelXY(w,wx,s,en,lvl)    | requires lvl between w[s] and w[en]
  555.     String w,wx
  556.     Variable/D s,en,lvl
  557.     
  558.     Variable/D p1,m
  559.     V_LevelX=FindLvlX($w,s,en,lvl)
  560.     if(exists(wx)==1)
  561.         p1=x2pnt($w,V_LevelX-(pnt2x($w,1)-pnt2x($w,0))/2)
  562.         m=(V_LevelX-pnt2x($w,p1))/(pnt2x($w,p1+1)-pnt2x($w,p1))
  563.         V_LevelX=$wx[p1] + m * ($wx[p1+1] -$wx[p1] )
  564.     endif
  565. End
  566.  
  567. |requires lvl between s and en, else returns NaN
  568. Function/D FindLvlX(w,s,en,lvl)
  569.     Wave/D w
  570.     Variable s,en |points
  571.     Variable/D lvl
  572.  
  573.     if(s==en)
  574.         return pnt2x(w,s)
  575.     endif
  576.     if((w[s]-lvl)*(w[en]-lvl) > 0)
  577.         return NaN
  578.     endif
  579.     Variable sp1,dir=sign(en-s)
  580.     Variable/D xpos,d,d1
  581.     | update s until s and s+dir bracket lvl
  582.     do
  583.         d=w[s]-lvl
  584.         if(d==0)
  585.             return pnt2x(w,s)
  586.         endif
  587.         sp1=s+dir
  588.         d1=w[sp1]-lvl
  589.         if(d1==0)
  590.             return pnt2x(w,sp1)
  591.         endif
  592.         if(d*d1<0) | on opposite sides of lvl
  593.             xpos=pnt2x(w,s)+(lvl-w[s])*(pnt2x(w,sp1)-pnt2x(w,s))/(w[sp1]-w[s])
  594.             break;
  595.         endif
  596.         s+=dir
  597.     while(s!=en) | same side
  598.     return xpos
  599. End
  600.  
  601. Proc AreaKind(kind,w,wx,s,en)    | outputs V_Area,S_Text
  602.     Variable kind,s,en |type, start,end points
  603.     String w,wx
  604.     Variable wlx=pnt2x($w,s),wrx=pnt2x($w,en)
  605.     if(kind==1)
  606.         S_Text="count"
  607.         V_Area=(en-s+1)*mean($w,wlx,wrx)
  608.     else
  609.     if(kind==2)
  610.         S_Text="rectangular area"
  611.         if(exists(wx)==1)
  612.             V_Area=AreaXYRPt($wx,$w,s,en)
  613.         else
  614.             V_Area=(wrx-wlx)*mean($w,wlx,wrx)
  615.         endif
  616.     else
  617.         S_Text="trapezoidal area"
  618.         if(exists(wx)==1)
  619.             V_Area=AreaXYTPt($wx,$w,s,en)
  620.         else
  621.             V_Area=area($w,wlx,wrx)
  622.         endif
  623.     endif
  624.     endif
  625. End
  626.  
  627. | returns number of knots in kn wave
  628. Function KnotsXYToBase(w,wx,wb,kn)
  629.     Wave w,wx,wb,kn
  630.     Variable knots=1
  631.     Variable p1=PrevNonNaN(kn,numpnts(w)) | right-most
  632.     Variable y1=w[p1]
  633.     Variable x1=wx[p1]
  634.     Variable p2,y2,x2,pt,m
  635.     do
  636.         p2=p1
  637.         p1=PrevNonNaN(kn,p2)
  638.         if(p1<0)
  639.             break;
  640.         endif
  641.         y2=y1
  642.         x2=x1
  643.         y1=w[p1]
  644.         x1=wx[p1]
  645.         pt=p1
  646.         if(x2==x1)
  647.             m=0
  648.         else
  649.             m=(y2-y1)/(x2-x1)
  650.         endif
  651.         do
  652.             wb[pt]=y1+m*(wx[pt]-x1)
  653.             pt+=1
  654.         while(pt<=p2)
  655.         knots+=1
  656.     while(1)    
  657.     return knots
  658. End
  659.  
  660. | returns number of knots in kn wave
  661. Function KnotsYToBase(w,wb,kn)
  662.     Wave w,wb,kn
  663.     Variable knots=1
  664.     Variable p1=PrevNonNaN(kn,numpnts(w)) | right-most
  665.     Variable y1=w[p1]
  666.     Variable x1=pnt2x(w,p1)
  667.     Variable p2,y2,x2,pt,m
  668.     do
  669.         p2=p1
  670.         p1=PrevNonNaN(kn,p2)
  671.         if(p1<0)
  672.             break;
  673.         endif
  674.         y2=y1
  675.         x2=x1
  676.         y1=w[p1]
  677.         x1=pnt2x(w,p1)
  678.         pt=p1
  679.         m=(y2-y1)/(x2-x1)
  680.         do
  681.             wb[pt]=y1+m*(pnt2x(w,pt)-x1)
  682.             pt+=1
  683.         while(pt<=p2)
  684.         knots+=1
  685.     while(1)    
  686.     return knots
  687. End
  688.  
  689.  
  690. Function PrevNonNaN(w,pt) | returns point number of previous nonNaN value in w
  691.     Wave w
  692.     Variable pt
  693.     do
  694.         pt-=1
  695.     while((pt>=0) %& (numtype(w[pt])==2))
  696.     
  697.     return pt | returns -1 if last
  698. end
  699.  
  700. Function NextNonNaN(w,pt) | returns point number of previous nonNaN value in w
  701.     Wave w
  702.     Variable pt
  703.     Variable lim=numpnts(w)
  704.     do
  705.         pt+=1
  706.     while((pt<lim) %& (numtype(w[pt])==2))
  707.     
  708.     return pt | returns numpnts(w) if last
  709. end
  710.  
  711. Proc XBoxToPBox(bx,w,wx,s,en)
  712.     Variable bx,s,en
  713.     String w,wx
  714.  
  715.     bx=abs(bx)
  716.     if(exists(wx)==1)
  717.         bx *= sign($wx[en]-$wx[s])
  718.         FindPX(w,wx,bx+$wx[s]);bx=V_peakP
  719.         FindPX(w,wx,$wx[s]);bx=abs(V_peakP-bx)
  720.     else
  721.         bx*=sign(rightx($w)-leftx($w))
  722.         bx=x2pnt($w,bx+leftx($w))
  723.     endif
  724.     Variable/G g_pbox=max(1,bx)
  725. End
  726.  
  727. Function/D MyBoxSmooth(w,pp,box)
  728.     Wave/D w
  729.     Variable pp,box | box should be odd
  730.     
  731.     Variable/D result=0
  732.     Variable  hp,top=numpnts(w)-1,terms
  733.     box = 2*trunc(box/2)+1
  734.     pp-=trunc(box/2)
  735.     hp=limit(0,pp+box-1,top)
  736.     pp=limit(0,pp,top)
  737.     terms=hp-pp+1
  738.     do
  739.         result+=w[pp]
  740.         pp+=1
  741.     while (pp<=hp)
  742.     return result/terms
  743. End
  744.  
  745.  
  746. Function AreaXYTPt(wx,wy,pa,pb)
  747.     Wave wx,wy    | monotonic wx!
  748.     Variable pa,pb    |pa<pb
  749.  
  750.     Variable a=0
  751.     do
  752.         a+=(wy[pa]+wy[pa+1])*(wx[pa+1]-wx[pa])/2
  753.         pa+=1
  754.     while(pa<pb)
  755.     return a
  756. End
  757.  
  758. Function AreaXYRPt(wx,wy,pa,pb)
  759.     Wave wx,wy    | monotonic wx!
  760.     Variable pa,pb    | pa<pb
  761.  
  762.     Variable a=0
  763.     do
  764.         a+=wy[pa]*(wx[pa+1]-wx[pa])
  765.         pa+=1
  766.     while(pa<pb)
  767.     return a
  768. End
  769.  
  770. | W_PM[0]=n, number of functions to fit
  771. | W_PM[1]=start p of peak fit
  772. | W_PM[2]=end p of peak fit
  773. | W_PM[3]=tp, type of baseline
  774. | W_PM[4]=nw, baseline terms, could be 0
  775. | W_PM[5]=tp, type of 1st peak function
  776. |...
  777. | w[0]=NaN
  778. | w[1]=K0
  779. | w[2]=1st  coefficient...
  780. Function/D PolyMorph(w,xx)
  781.     Wave/D w
  782.     Variable/D xx
  783.     
  784.     Variable/D tp,nw, bs=3,n=W_PM[0],st=2,ii,s,r=w[1] | K0
  785.      if(n==0)
  786.         return r
  787.     endif
  788.     do    
  789.         tp=W_PM[bs]
  790.         nw=W_PM[bs+1]
  791.         if(nw>0)
  792.             if(tp<4) | line, poly 3 - poly 5 K1*x+K2*x^2...
  793.                 ii=nw-1
  794.                 s=0
  795.                 do
  796.                     s=w[ii+st]+xx*s
  797.                     ii-=1
  798.                 while (ii >= 0)
  799.                 r+= s*xx
  800.             else
  801.                 if(tp==8) | gauss
  802.                     r+=w[st]*exp(-((xx-w[st+1])/w[st+2])^2)
  803.                 else
  804.                 if(tp==7) | lor
  805.                     r+=w[st]/((xx-w[st+1])^2+w[st+2])
  806.                 else
  807.                 if(tp==6) | exp
  808.                     r+= w[st]*exp(-w[st+1]*xx)
  809.                 else
  810.                 if(tp==5)    | dblexp
  811.                     r+= w[st]*exp(-w[st+1]*xx)*exp(-w[st+2]*xx)
  812.                 else    | 4, sin
  813.                     r+= w[st]*sin(w[st+1]*xx+w[st+2])
  814.                 endif
  815.                 endif
  816.                 endif
  817.                 endif
  818.             endif
  819.             st+=nw    
  820.         endif
  821.         bs+=2
  822.         n-=1
  823.     while (n>0)
  824.     return r
  825. End
  826.  
  827.  
  828. Proc FindPX(w,wx,xx)
  829.     String w,wx    | wx monotonic
  830.     Variable xx
  831.  
  832.     if(exists(wx)==1)
  833.         FindLevel/Q $wx, xx 
  834.         V_peakP = x2pnt($wx,V_levelX)
  835.         if(V_Flag==1)
  836.             if((xx<$wx[1])%&($wx[0]<$wx[1]))
  837.                 V_peakP=0
  838.             else
  839.                 V_peakP=numpnts($w)-1
  840.             endif        
  841.         endif
  842.         V_peakX = pnt2x($w,V_peakP)
  843.     else
  844.         V_peakP = x2pnt($w,xx)
  845.         V_peakX = xx
  846.     endif
  847. End
  848.  
  849. Proc ChkLen(w1,w2)
  850.     String w1,w2    | Wave names
  851.     if((exists(w1)==1) %& (exists(w2)==1))
  852.         if(numpnts($w1) != numpnts($w2))
  853.             Abort w1+" and "+w2+" have different number of points!"
  854.         endif
  855.     endif
  856. End
  857.  
  858. Proc Mk(w,n)
  859.     String w
  860.     Variable n
  861.     if(exists(w)==1)
  862.         Redimension/D/N=(n) $w;$w=NaN
  863.     else
  864.         Make/D/N=(n) $w=NaN
  865.     endif    
  866. End
  867.  
  868. Proc Dup(t,w)
  869.     String t,w
  870.     String n=""
  871.     if (exists(w)==1)
  872.         n=note($w) | keep w's note
  873.     endif
  874.     Duplicate/O $t,$w
  875.     Note/K $w;Note $w,n
  876. End    
  877.  
  878. Proc SameLen(t,w)
  879.     String t,w
  880.     if(exists(w)!=1)
  881.         Duplicate/O $t,$w
  882.     else
  883.         Redimension/D/N=(numpnts($t)) $w;CopyScales $t,$w
  884.     endif
  885. End    
  886.  
  887. Proc NewWv(srcw,sfx)
  888.     String srcw,sfx
  889.  
  890.     Variable ii=1
  891.     S_Wave=srcw[0,17-strlen(sfx)]+sfx
  892.     if(exists(S_Wave)==1)
  893.         String shw=S_Wave
  894.         do
  895.             S_Wave=srcw[0,14-strlen(sfx)]+num2istr(ii)+sfx
  896.             ii+=1
  897.         while ((exists(S_Wave)==1)*(ii<999))
  898.     endif
  899.     if(exists(srcw)==1)
  900.         Duplicate/O $srcw,$S_Wave
  901.     else
  902.         Make/D/O $S_Wave
  903.     endif
  904.     Print "Created wave "+S_Wave
  905. End
  906.  
  907. Proc RmWv(w)
  908.     String w
  909.  
  910.     CheckDisplayed/W=$WinName(0,1) $w
  911.     if(V_Flag==1)
  912.         DoWindow/F $WinName(0,1);Remove $w
  913.     endif
  914. End
  915.  
  916. Proc KillWv(w)
  917.     String w
  918.     if(exists(w)==1)
  919.         CheckDisplayed/A $w
  920.         if(V_Flag==0)
  921.             KillWaves $w
  922.         endif
  923.     endif
  924. End
  925.  
  926. Proc AppCrs(w)
  927.     String w
  928.     ShowInfo;Cursor/P A,$w,0;Cursor/P B,$w,numpnts($w)-1
  929. End
  930. Proc NoCrs()
  931.     HideInfo;Cursor/K A;Cursor/K B
  932. End
  933. Proc AppWv(w,wx,ax)
  934.     String w,wx,ax
  935.     
  936.     String wn=WinName(0,1)
  937.     if(strlen(wn)==0)
  938.         ax = "Display $w"
  939.     else
  940.         DoWindow/F $wn
  941.         ax="Append"+ax+" $w "
  942.     endif
  943.     if(exists(w)==1)
  944.         CheckDisplayed $w
  945.         if(V_Flag==0)
  946.             if(exists(wx)==1)
  947.                 ax+=" vs $wx"
  948.             endif
  949.             Execute ax
  950.             ColorWaves()
  951.         endif
  952.     endif
  953. End
  954. Proc SBs(s)
  955.     String s
  956.     if((exists(s)==1)%|(cmpstr(s+";",none)==0))
  957.         g_b=s;rb_b=s;fpks_b=s;afp_b=s
  958.     else
  959.     if(strsearch(S_funcs,s,0)>=0)
  960.         fpks_b=s
  961.     else
  962.     if(cmpstr(s,"_From Fit Peaks_")==0)
  963.         rb_b=s;afp_b=s
  964.     else
  965.     if(cmpstr(s,"_From Fit Baseline_")==0)
  966.         fpks_b=s
  967.     endif
  968.     endif
  969.     endif
  970.     endif
  971. End
  972.  
  973. Macro ShowOnlyDataAndBase()
  974.     Silent 1;PauseUpdate
  975.     String keep=g_w+";"+g_wx+";"
  976.     keep+=g_b+";"+WaveList("*X",";","WIN:"+WinName(0,1))
  977.     keep+=g_keep+";";g_keep=""
  978.     String w,wn=WinName(0,1),ws=WaveList("*",";","WIN:"+wn)
  979.     Variable pos,flg=0
  980.     if(strlen(wn)!=0)
  981.         do
  982.             pos=strsearch(ws,";",0)
  983.             if(pos<0)
  984.                 break
  985.             endif
  986.             w=ws[0,pos-1]
  987.             ws[0,pos]=""
  988.             if(strsearch(keep,w+";",0)<0)
  989.                 RmWv(w);flg=1
  990.             endif
  991.         while (1)
  992.         if(flg)
  993.             ColorWaves()
  994.         endif
  995.     endif
  996. End
  997.  
  998. Macro ColorWaves() : GraphStyle
  999.     Silent 1;PauseUpdate
  1000.     Modify/Z rgb[0]=(65535,0,0)
  1001.     Modify/Z rgb[1]=(0,0,65535),rgb[2]=(3009,65535,1882),rgb[3]=(0,0,0)
  1002.     Modify/Z rgb[4]=(0,0,65535),rgb[5]=(63953,3661,65535)
  1003.     Modify/Z rgb[6]=(37510,1,1),rgb[7]=(27232,40528,22540),rgb[8]=(1531,1314,28456)
  1004.     Legend/C/N=default ""
  1005. End
  1006.  
  1007.  
  1008. Proc CheckTwoCursors(w)
  1009.     String w    | the y wave
  1010.     Variable/G V_start,V_theEnd,V_startX,V_theEndX | point indices, x vals
  1011.     String wa=CsrWave(A),wb=CsrWave(B),t="cursors on target graph "
  1012.     if((strlen(wa)==0)%|(strlen(wb)==0))
  1013.         Abort "expecting two "+t
  1014.     endif
  1015.     V_start=pcsr(A);V_theEnd=pcsr(B)
  1016.     V_startX=hcsr(A);V_theEndX=hcsr(B)
  1017.     Variable x0=pnt2x($w,0),x1=pnt2x($w,1),ok
  1018.     ok = (numpnts($wa)==numpnts($w)) %& (numpnts($wb)==numpnts($w)) 
  1019.     ok = ok %& (x0==pnt2x($wa,0)) %& (x1==pnt2x($wa,1))
  1020.     ok = ok %& (x0==pnt2x($wb,0)) %& (x1==pnt2x($wb,1)) 
  1021.     if(ok)
  1022.         if(V_start>V_theEnd)
  1023.             x0=V_Start
  1024.             V_start=V_theEnd
  1025.             V_theEnd=x0
  1026.             x0=V_startX
  1027.             V_startX=V_theEndX
  1028.             V_theEndX=x0
  1029.         endif
  1030. |        if(V_start==V_theEnd)
  1031. |            Abort t+"are too close together!"
  1032. |        endif
  1033.     else
  1034.         Abort t+"must be on wave "+w+" or one like it"
  1035.     endif
  1036. End
  1037.  
  1038. Window PeakReportTable() : Table
  1039.     PauseUpdate; Silent 1        | building window...
  1040.     Edit /W=(14,46,545,227) W_PkCenters.y,W_PkAmps.y,W_PkFWHMs.y as "PeakReportTable"
  1041.     Append W_PkAreas.y,W_PkX1.y,W_PkX2.y,W_PkXwL.y
  1042.     Modify size=9
  1043.     Modify width(Point)=38
  1044.     Modify width(W_PkCenters.y)=46
  1045.     Modify width(W_PkAmps.y)=72
  1046.     Modify width(W_PkFWHMs.y)=78
  1047.     Modify width(W_PkAreas.y)=72
  1048.     Modify width(W_PkX1.y)=50
  1049.     Modify width(W_PkX2.y)=52
  1050.     Modify width(W_PkXwL.y)=94
  1051. EndMacro
  1052.  
  1053. Window PeakReportLayout() : Layout
  1054.     PauseUpdate; Silent 1        | building window...
  1055.     Layout /W=(18,45,536,265) PeakReportTable(49,190,553,310)/O=2 as "Peak Report Layout"
  1056.     Textbox /N=title/A=LT/X=18.5455/Y=2.06044 "\JCYour Report Title Here"
  1057.     Modify left(title)=133,top(title)=47,width(title)=105,height(title)=4,frame(title)=4
  1058.     Textbox /N=info/A=LT/X=63.8182/Y=0.961538 "7 peaks\r(trapezoidal area)"
  1059.     Modify left(info)=382,top(info)=39,width(info)=83,height(info)=18,frame(info)=1
  1060.     Append PeakFitGraph(43,64,555,183)/O=1
  1061. EndMacro
  1062.  
  1063.  
  1064.  
  1065.  
  1066.  
  1067. Window PeakEstimates() : Table
  1068.     Edit /W=(35,87,576,268) W_EstCentersP.y,W_EstCentersX.y,W_EstAmpsY.y as "Peak Estimates"
  1069.     Append W_EstWidthsX.y,W_EstEdgesP.y
  1070.     Modify width=98,width(Point)=64
  1071. EndMacro
  1072.  
  1073. Window AreaReportTable() : Table
  1074.     Edit /W=(59,222,420,347) W_AreaX1.y,W_AreaX2.y,W_AreaNoBase.y as "Areas Between Cursors"
  1075.     Modify width(Point)=64
  1076.     Modify width(W_AreaNoBase.y)=112
  1077. EndMacro
  1078.  
  1079.